home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / ode-initval / rk8pd.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  10.7 KB  |  434 lines
  1. /* ode-initval/rk8pd.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Runge-Kutta 8(9), Prince-Dormand */
  21.  
  22. /* Author:  G. Jungman
  23.  */
  24. #include <config.h>
  25. #include <stdlib.h>
  26. #include <string.h>
  27. #include <gsl/gsl_errno.h>
  28. #include <gsl/gsl_odeiv.h>
  29.  
  30. #include "odeiv_util.h"
  31.  
  32. /* Prince-Dormand constants */
  33.  
  34. static const double Abar[] = {
  35.   14005451.0 / 335480064.0,
  36.   0.0,
  37.   0.0,
  38.   0.0,
  39.   0.0,
  40.   -59238493.0 / 1068277825.0,
  41.   181606767.0 / 758867731.0,
  42.   561292985.0 / 797845732.0,
  43.   -1041891430.0 / 1371343529.0,
  44.   760417239.0 / 1151165299.0,
  45.   118820643.0 / 751138087.0,
  46.   -528747749.0 / 2220607170.0,
  47.   1.0 / 4.0
  48. };
  49.  
  50. static const double A[] = {
  51.   13451932.0 / 455176623.0,
  52.   0.0,
  53.   0.0,
  54.   0.0,
  55.   0.0,
  56.   -808719846.0 / 976000145.0,
  57.   1757004468.0 / 5645159321.0,
  58.   656045339.0 / 265891186.0,
  59.   -3867574721.0 / 1518517206.0,
  60.   465885868.0 / 322736535.0,
  61.   53011238.0 / 667516719.0,
  62.   2.0 / 45.0
  63. };
  64.  
  65. static const double ah[] = {
  66.   1.0 / 18.0,
  67.   1.0 / 12.0,
  68.   1.0 / 8.0,
  69.   5.0 / 16.0,
  70.   3.0 / 8.0,
  71.   59.0 / 400.0,
  72.   93.0 / 200.0,
  73.   5490023248.0 / 9719169821.0,
  74.   13.0 / 20.0,
  75.   1201146811.0 / 1299019798.0
  76. };
  77.  
  78. static const double b21 = 1.0 / 18.0;
  79. static const double b3[] = { 1.0 / 48.0, 1.0 / 16.0 };
  80. static const double b4[] = { 1.0 / 32.0, 0.0, 3.0 / 32.0 };
  81. static const double b5[] = { 5.0 / 16.0, 0.0, -75.0 / 64.0, 75.0 / 64.0 };
  82. static const double b6[] = { 3.0 / 80.0, 0.0, 0.0, 3.0 / 16.0, 3.0 / 20.0 };
  83. static const double b7[] = {
  84.   29443841.0 / 614563906.0,
  85.   0.0,
  86.   0.0,
  87.   77736538.0 / 692538347.0,
  88.   -28693883.0 / 1125000000.0,
  89.   23124283.0 / 1800000000.0
  90. };
  91. static const double b8[] = {
  92.   16016141.0 / 946692911.0,
  93.   0.0,
  94.   0.0,
  95.   61564180.0 / 158732637.0,
  96.   22789713.0 / 633445777.0,
  97.   545815736.0 / 2771057229.0,
  98.   -180193667.0 / 1043307555.0
  99. };
  100. static const double b9[] = {
  101.   39632708.0 / 573591083.0,
  102.   0.0,
  103.   0.0,
  104.   -433636366.0 / 683701615.0,
  105.   -421739975.0 / 2616292301.0,
  106.   100302831.0 / 723423059.0,
  107.   790204164.0 / 839813087.0,
  108.   800635310.0 / 3783071287.0
  109. };
  110. static const double b10[] = {
  111.   246121993.0 / 1340847787.0,
  112.   0.0,
  113.   0.0,
  114.   -37695042795.0 / 15268766246.0,
  115.   -309121744.0 / 1061227803.0,
  116.   -12992083.0 / 490766935.0,
  117.   6005943493.0 / 2108947869.0,
  118.   393006217.0 / 1396673457.0,
  119.   123872331.0 / 1001029789.0
  120. };
  121. static const double b11[] = {
  122.   -1028468189.0 / 846180014.0,
  123.   0.0,
  124.   0.0,
  125.   8478235783.0 / 508512852.0,
  126.   1311729495.0 / 1432422823.0,
  127.   -10304129995.0 / 1701304382.0,
  128.   -48777925059.0 / 3047939560.0,
  129.   15336726248.0 / 1032824649.0,
  130.   -45442868181.0 / 3398467696.0,
  131.   3065993473.0 / 597172653.0
  132. };
  133. static const double b12[] = {
  134.   185892177.0 / 718116043.0,
  135.   0.0,
  136.   0.0,
  137.   -3185094517.0 / 667107341.0,
  138.   -477755414.0 / 1098053517.0,
  139.   -703635378.0 / 230739211.0,
  140.   5731566787.0 / 1027545527.0,
  141.   5232866602.0 / 850066563.0,
  142.   -4093664535.0 / 808688257.0,
  143.   3962137247.0 / 1805957418.0,
  144.   65686358.0 / 487910083.0
  145. };
  146. static const double b13[] = {
  147.   403863854.0 / 491063109.0,
  148.   0.0,
  149.   0.0,
  150.   -5068492393.0 / 434740067.0,
  151.   -411421997.0 / 543043805.0,
  152.   652783627.0 / 914296604.0,
  153.   11173962825.0 / 925320556.0,
  154.   -13158990841.0 / 6184727034.0,
  155.   3936647629.0 / 1978049680.0,
  156.   -160528059.0 / 685178525.0,
  157.   248638103.0 / 1413531060.0,
  158.   0.0
  159. };
  160.  
  161. typedef struct
  162. {
  163.   double *k[13];
  164.   double *ytmp;
  165. }
  166. rk8pd_state_t;
  167.  
  168. static void *
  169. rk8pd_alloc (size_t dim)
  170. {
  171.   rk8pd_state_t *state = (rk8pd_state_t *) malloc (sizeof (rk8pd_state_t));
  172.   int i, j;
  173.  
  174.   if (state == 0)
  175.     {
  176.       GSL_ERROR_NULL ("failed to allocate space for rk8pd_state", GSL_ENOMEM);
  177.     }
  178.  
  179.   state->ytmp = (double *) malloc (dim * sizeof (double));
  180.  
  181.   if (state->ytmp == 0)
  182.     {
  183.       free (state);
  184.       GSL_ERROR_NULL ("failed to allocate space for ytmp", GSL_ENOMEM);
  185.     }
  186.  
  187.   for (i = 0; i < 13; i++)
  188.     {
  189.       state->k[i] = (double *) malloc (dim * sizeof (double));
  190.  
  191.       if (state->k[i] == 0)
  192.     {
  193.       for (j = 0; j < i; j++)
  194.         {
  195.           free (state->k[j]);
  196.         }
  197.       free (state->ytmp);
  198.       free (state);
  199.       GSL_ERROR_NULL ("failed to allocate space for k's", GSL_ENOMEM);
  200.     }
  201.     }
  202.  
  203.   return state;
  204. }
  205.  
  206.  
  207. static int
  208. rk8pd_apply (void *vstate,
  209.          size_t dim,
  210.          double t,
  211.          double h,
  212.          double y[],
  213.          double yerr[],
  214.          const double dydt_in[],
  215.          double dydt_out[], const gsl_odeiv_system * sys)
  216. {
  217.   rk8pd_state_t *state = (rk8pd_state_t *) vstate;
  218.  
  219.   size_t i;
  220.   int status = 0;
  221.  
  222.   double *const ytmp = state->ytmp;
  223.   /* Note that k1 is stored in state->k[0] due to zero-based indexing */
  224.   double *const k1 = state->k[0];
  225.   double *const k2 = state->k[1];
  226.   double *const k3 = state->k[2];
  227.   double *const k4 = state->k[3];
  228.   double *const k5 = state->k[4];
  229.   double *const k6 = state->k[5];
  230.   double *const k7 = state->k[6];
  231.   double *const k8 = state->k[7];
  232.   double *const k9 = state->k[8];
  233.   double *const k10 = state->k[9];
  234.   double *const k11 = state->k[10];
  235.   double *const k12 = state->k[11];
  236.   double *const k13 = state->k[12];
  237.  
  238.   /* k1 step */
  239.   if (dydt_in != NULL)
  240.     {
  241.       DBL_MEMCPY (k1, dydt_in, dim);
  242.     }
  243.   else
  244.     {
  245.       int s = GSL_ODEIV_FN_EVAL (sys, t, y, k1);
  246.       GSL_STATUS_UPDATE (&status, s);
  247.     }
  248.  
  249.   for (i = 0; i < dim; i++)
  250.     ytmp[i] = y[i] + b21 * h * k1[i];
  251.  
  252.   /* k2 step */
  253.   {
  254.     int s = GSL_ODEIV_FN_EVAL (sys, t + ah[0] * h, ytmp, k2);
  255.     GSL_STATUS_UPDATE (&status, s);
  256.   }
  257.   for (i = 0; i < dim; i++)
  258.     ytmp[i] = y[i] + h * (b3[0] * k1[i] + b3[1] * k2[i]);
  259.  
  260.   /* k3 step */
  261.   {
  262.     int s = GSL_ODEIV_FN_EVAL (sys, t + ah[1] * h, ytmp, k3);
  263.     GSL_STATUS_UPDATE (&status, s);
  264.   }
  265.   for (i = 0; i < dim; i++)
  266.     ytmp[i] = y[i] + h * (b4[0] * k1[i] + b4[2] * k3[i]);
  267.  
  268.   /* k4 step */
  269.   {
  270.     int s = GSL_ODEIV_FN_EVAL (sys, t + ah[2] * h, ytmp, k4);
  271.     GSL_STATUS_UPDATE (&status, s);
  272.   }
  273.   for (i = 0; i < dim; i++)
  274.     ytmp[i] = y[i] + h * (b5[0] * k1[i] + b5[2] * k3[i] + b5[3] * k4[i]);
  275.  
  276.   /* k5 step */
  277.   {
  278.     int s = GSL_ODEIV_FN_EVAL (sys, t + ah[3] * h, ytmp, k5);
  279.     GSL_STATUS_UPDATE (&status, s);
  280.   }
  281.   for (i = 0; i < dim; i++)
  282.     ytmp[i] = y[i] + h * (b6[0] * k1[i] + b6[3] * k4[i] + b6[4] * k5[i]);
  283.  
  284.   /* k6 step */
  285.   {
  286.     int s = GSL_ODEIV_FN_EVAL (sys, t + ah[4] * h, ytmp, k6);
  287.     GSL_STATUS_UPDATE (&status, s);
  288.   }
  289.   for (i = 0; i < dim; i++)
  290.     ytmp[i] =
  291.       y[i] + h * (b7[0] * k1[i] + b7[3] * k4[i] + b7[4] * k5[i] +
  292.           b7[5] * k6[i]);
  293.  
  294.   /* k7 step */
  295.   {
  296.     int s = GSL_ODEIV_FN_EVAL (sys, t + ah[5] * h, ytmp, k7);
  297.     GSL_STATUS_UPDATE (&status, s);
  298.   }
  299.   for (i = 0; i < dim; i++)
  300.     ytmp[i] =
  301.       y[i] + h * (b8[0] * k1[i] + b8[3] * k4[i] + b8[4] * k5[i] +
  302.           b8[5] * k6[i] + b8[6] * k7[i]);
  303.  
  304.   /* k8 step */
  305.   {
  306.     int s = GSL_ODEIV_FN_EVAL (sys, t + ah[6] * h, ytmp, k8);
  307.     GSL_STATUS_UPDATE (&status, s);
  308.   }
  309.   for (i = 0; i < dim; i++)
  310.     ytmp[i] =
  311.       y[i] + h * (b9[0] * k1[i] + b9[3] * k4[i] + b9[4] * k5[i] +
  312.           b9[5] * k6[i] + b9[6] * k7[i] + b9[7] * k8[i]);
  313.  
  314.   /* k9 step */
  315.   {
  316.     int s = GSL_ODEIV_FN_EVAL (sys, t + ah[7] * h, ytmp, k9);
  317.     GSL_STATUS_UPDATE (&status, s);
  318.   }
  319.   for (i = 0; i < dim; i++)
  320.     ytmp[i] =
  321.       y[i] + h * (b10[0] * k1[i] + b10[3] * k4[i] + b10[4] * k5[i] +
  322.           b10[5] * k6[i] + b10[6] * k7[i] + b10[7] * k8[i] +
  323.           b10[8] * k9[i]);
  324.  
  325.   /* k10 step */
  326.   {
  327.     int s = GSL_ODEIV_FN_EVAL (sys, t + ah[8] * h, ytmp, k10);
  328.     GSL_STATUS_UPDATE (&status, s);
  329.   }
  330.   for (i = 0; i < dim; i++)
  331.     ytmp[i] =
  332.       y[i] + h * (b11[0] * k1[i] + b11[3] * k4[i] + b11[4] * k5[i] +
  333.           b11[5] * k6[i] + b11[6] * k7[i] + b11[7] * k8[i] +
  334.           b11[8] * k9[i] + b11[9] * k10[i]);
  335.  
  336.   /* k11 step */
  337.   {
  338.     int s = GSL_ODEIV_FN_EVAL (sys, t + ah[9] * h, ytmp, k11);
  339.     GSL_STATUS_UPDATE (&status, s);
  340.   }
  341.   for (i = 0; i < dim; i++)
  342.     ytmp[i] =
  343.       y[i] + h * (b12[0] * k1[i] + b12[3] * k4[i] + b12[4] * k5[i] +
  344.           b12[5] * k6[i] + b12[6] * k7[i] + b12[7] * k8[i] +
  345.           b12[8] * k9[i] + b12[9] * k10[i] + b12[10] * k11[i]);
  346.  
  347.   /* k12 step */
  348.   {
  349.     int s = GSL_ODEIV_FN_EVAL (sys, t + h, ytmp, k12);
  350.     GSL_STATUS_UPDATE (&status, s);
  351.   }
  352.   for (i = 0; i < dim; i++)
  353.     ytmp[i] =
  354.       y[i] + h * (b13[0] * k1[i] + b13[3] * k4[i] + b13[4] * k5[i] +
  355.           b13[5] * k6[i] + b13[6] * k7[i] + b13[7] * k8[i] +
  356.           b13[8] * k9[i] + b13[9] * k10[i] + b13[10] * k11[i] +
  357.           b13[11] * k12[i]);
  358.  
  359.   /* k13 step */
  360.   {
  361.     int s = GSL_ODEIV_FN_EVAL (sys, t + h, ytmp, k13);
  362.     GSL_STATUS_UPDATE (&status, s);
  363.   }
  364.  
  365.   /* final sum and error estimate */
  366.   for (i = 0; i < dim; i++)
  367.     {
  368.       const double ksum8 =
  369.     Abar[0] * k1[i] + Abar[5] * k6[i] + Abar[6] * k7[i] +
  370.     Abar[7] * k8[i] + Abar[8] * k9[i] + Abar[9] * k10[i] +
  371.     Abar[10] * k11[i] + Abar[11] * k12[i] + Abar[12] * k13[i];
  372.       const double ksum7 =
  373.     A[0] * k1[i] + A[5] * k6[i] + A[6] * k7[i] + A[7] * k8[i] +
  374.     A[8] * k9[i] + A[9] * k10[i] + A[10] * k11[i] + A[11] * k12[i];
  375.       y[i] += h * ksum8;
  376.       yerr[i] = h * (ksum7 - ksum8);
  377.       if (dydt_out != NULL)
  378.     dydt_out[i] = ksum8;
  379.     }
  380.  
  381.   return status;
  382. }
  383.  
  384. static int
  385. rk8pd_reset (void *vstate, size_t dim)
  386. {
  387.   rk8pd_state_t *state = (rk8pd_state_t *) vstate;
  388.  
  389.   int i;
  390.  
  391.   for (i = 0; i < 13; i++)
  392.     {
  393.       DBL_ZERO_MEMSET (state->k[i], dim);
  394.     }
  395.  
  396.   DBL_ZERO_MEMSET (state->ytmp, dim);
  397.  
  398.   return GSL_SUCCESS;
  399. }
  400.  
  401. static unsigned int
  402. rk8pd_order (void *vstate)
  403. {
  404.   rk8pd_state_t *state = (rk8pd_state_t *) vstate;
  405.   state = 0; /* prevent warnings about unused parameters */
  406.   return 8;
  407. }
  408.  
  409. static void
  410. rk8pd_free (void *vstate)
  411. {
  412.   rk8pd_state_t *state = (rk8pd_state_t *) vstate;
  413.   int i;
  414.  
  415.   for (i = 0; i < 13; i++)
  416.     {
  417.       free (state->k[i]);
  418.     }
  419.   free (state->ytmp);
  420.   free (state);
  421. }
  422.  
  423. static const gsl_odeiv_step_type rk8pd_type = { "rk8pd",    /* name */
  424.   1,                /* can use dydt_in */
  425.   0,                /* gives exact dydt_out */
  426.   &rk8pd_alloc,
  427.   &rk8pd_apply,
  428.   &rk8pd_reset,
  429.   &rk8pd_order,
  430.   &rk8pd_free
  431. };
  432.  
  433. const gsl_odeiv_step_type *gsl_odeiv_step_rk8pd = &rk8pd_type;
  434.